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Abstract: Many drug delivery systems suffer from undesirable interactions with the host im- 
mune system. It has been experimentally established that covalent attachment (irreversible 
adsorption) of suitable macromolecules to the surface of the drug carrier can reduce such unde- 
sirable interactions. A fundamental understanding of the adsorption process is still lacking. In 
this paper, the classical random irreversible adsorption model is generalized to capture certain 
essential processes involved in pharmacological applications, allowing for macromolecules of 
different sizes, partial overlapping of the tails of macromolecules, and the influence of reactions 
with the solvent on the adsorption process. Working in one dimension, an integro-differential 
evolution equation for the adsorption process is derived and the asymptotic behaviour of the 
surface area covered and the number of molecules attached to the surface is studied. Fi- 
nally, equation-free dynamic renormalization tools are applied to study the asymptotically 
self-similar behaviour of the adsorption statistics. 



1 Introduction 

Random sequential adsorption (RSA) is a classical model for various physical, chemical 
or biological problems [Sj. In the simplest form, RSA processes can be formulated as se- 
quential addition to a structure of objects that cannot overlap, and once inserted, cannot 
move or leave the structure jT| . In this paper, we present a pharmacological example in 
which application of the RSA model can provide meaningful qualitative insights. Moti- 
vated by pharmacological applications, we present a slight generalization of the classical 
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RSA model to enable us to study the effects of polydispersity and partial overlap of 
adsorbing macromolecules on the surface of a virus. We also study the dependence of 
the adsorption process on interactions (reactions) of the adsorbing macromolecules with 
the solvent. 

The paper is organized as follows. In Section |21 we will introduce the motivating 
pharmacological example and our questions of interest. In Section H3 we will introduce 
the generalized random sequential adsorption (gRSA) model suitable for capturing es- 
sential features of the pharmacological problem from Section which we formulate in 
one dimension. The analytical results for this model are presented in Section HJ We 
derive the governing integro-differential equation for the evolution of gaps between poly- 
mers, and compute the asymptotical properties of the quantities of interest, namely the 
number of macromolecules adsorbed and the total area (in one dimension, length) they 
cover. In Section El we apply equation-free methods to the computational study of the 
system. The main idea that underlies this equation-free computer-assisted analysis is the 
design and execution of appropriately-initialized short bursts of stochastic simulations; 
the results of these are processed to estimate coarse-grained quantities of interest - in 
this case the self-similarly evolving shape of the gap statistics in the problem. Finally, 
in Section we discuss the higher dimensional case and summarize the connections 
between the theory and the experimental data. 

2 Pharmacological background 

Many medical conditions such as cancer, heart disease and heritable disorders (hemophi- 
lia, cystic fibrosis, muscular dystrophy etc) have faulty, mutant genes as an underlying 
cause. Healthy, normal genes can be readily synthesized in the laboratory but introduc- 
ing them into diseased cells remains a challenge. Many research groups are studying 
viruses, such as adenovirus, as a means to introduce normal genes into diseased cells. 
For therapeutic use, the virus' own DNA is usually partially or completely replaced 
by the gene of interest. The most common adenovirus strain used for this purpose 
is adenovirus type 5 (Ad5), because it is easy to manipulate and is non-pathogenic 
in humans [Ti] . Ad5 has been used with great success to treat diseases in labora- 
tory animals but the results have not been replicated in humans. One of the greatest 
problems with using Ad5 in humans is the presence of neutralizing antibodies. In ad- 
dition, the viruses often infect non-target cells, particularly the liver, causing unwanted 
toxicity. Our laboratory is exploring the use of hydrophilic polymers such as poly[N- 
(2-hydroxypropyl)methacrylamide] (pHPMA) to coat virus particles and protect them 
from neutralising antibodies by steric shielding. This technique is very effective at pro- 
tecting the virus and permitting it to be retargeted to specific cell types El U2] • The 
polymer has multiple esters along its length that are used to bind to the amino side 
chain of lysine residues on the virus surface. In a coating reaction the polymers bind 
randomly to the virus surface until (a) all of the lysine residues are occupied or (b) 
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lysine residues are rendered inaccessible (obscured) by polymer chains. We know little 
about the orientation of polymers on the virus surface or how to optimise the coating 
reaction because there are no techniques to visualise the orientation of polymers on the 
virus surface. 

Representative questions one would like to answer are: How many polymer molecules 
will become attached to the viral surface by a given time? How large is the surface area 
covered by the polymer coat at that time? In this paper, the theoretical approach is 
chosen to address these questions for simplified models of the adsorption process. Since 
the adsorption process is driven by the diffusion of molecules to the surface of the virus, 
and since the adsorption is effectively irreversible, a suitable modification of the classical 
RSA might be applied to model the process. It is important to take into account poly- 
dispersity in the polymer size. Even if we prepare the polymer molecules with a specified 
target molecular weight, some relatively small molecules of the polymer will always be 
present, and they will diffuse faster than the larger molecules. A smaller molecule can 
reach the surface of the virus at a higher rate. We will consider the adsorption of poly- 
mers whose diameters are distributed according to a probability distribution function 
P(z). The solution is assumed to be well- mixed. We also define p(z) as the probability 
distribution function of the particles which can reach the surface in a single time step. 
An important modelling issue lies in a good choice of P(z) (e.g. it might be the Gamma 
distribution) and in a realistic relation between P(z) and p(z). In Sectional we simply 
specify p(z) (avoiding the above questions). The long term dynamics of the polydis- 
perse adsorption depends on the behaviour of p(z) close to zero. We will use different 
distributions p(z) given by (jH.ljl . 

Finally, the reactive groups of the polymer molecules can also react with the solvent 
before reaching the surface. This is the case for commonly used biocompatible pHPMA 
reactive polymers ^3] . If all reactive groups of a polymer are already hydrolyzed then 
the polymer cannot covalently bind to the surface of the virus. Hence, we have to 
consider that only a fraction r(t) of the polymers is still reactive at time t. Depending 
on the form of r(i), different polymer coats may be created. This issue will be discussed 
in more detail in Section [Ol 

Polymers are long flexible molecules 0]. The pHPMA polymer molecule has (one 
or more) reactive group (s) which can react with the primary amino groups on the viral 
surface. As a result, a polymer molecule becomes covalently (irreversibly) attached to 
the surface at a point. The rest of the polymer is not attached (unless another covalent 
bond is created) and it freely "wiggles" in the space above the viral surface. Having in 
mind that the "wiggling tail" does not perfectly shield the underlying surface, we gener- 
alize the classical RSA model to allow partial overlap of the adsorbing macromolecules, 
i.e. we allow some squeezing of the polymers. 

Let N(t) be the number of polymers attached to the surface at time t. Let A(t) be 
the total area of the surface covered by adsorbed polymers at time t. Since we cover 
the surface by polymers of different sizes, there is no obvious relation between A(t) and 
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N(t). However, both variables A(t) and N(t) are of practical interest as discussed below. 

We cover the surface of the virus by polymers to protect the surface from unwanted 
interactions. Hence, the number A(t) gives us the simple characterization of the area of 
the surface which is protected by the polymer coat. The unwanted interactions are not 
the only problem which one has to overcome in order to use viruses as a drug delivery 
system. Another important task is to retarget the virus to infect the cells of interest 
(e.g. cancer cells) via new receptors. Assuming that we put one "targeting" group per 
polymer molecule, the number of targeting molecules will be equal to N(t). 

If we considered the adsorption of the same-size nonoverlapping objects of the area 
a, then we would have A(t) = aN(t). In our case, the adsorbing molecules have different 
sizes. There is no obvious relation between A(t) and N(t) and both quantities are of 
interest. In the following sections, we will present theoretical approaches to compute 
the time evolution of A(t) and N(t). 

3 A simple generalization of random sequential ad- 
sorption 

Random sequential adsorption has been extensively studied during the last several 
decades The theoretical work is more mature in one dimension with information in 
higher dimensions mostly coming from numerical simulations jTUj. If we consider fixed 
size objects, then RSA usually starts from an empty surface and continues until the 
time when no further object can be placed, the so-called "jamming limit". If the objects 
to be covered have spherical symmetry ^3] then the coverage approaches the jamming 
limit as t~ x l d where t is time and d is the dimension. The asymptotic behaviour can be 
more complicated for objects of different shape [5]. 

As argued in Section [21 polydispersity is often present in real systems. If we allow ad- 
sorbing particles (in one spatial dimension) of arbitrarily small length, then the coverage 
approaches the full coverage as t — > oo. Relatively less is known for polydisperse adsorp- 
tion. One-dimensional analytical results can be found in jTU], where it is assumed that 
the attached polymers prevent binding of other polymer molecules that would overlap 
with them. In reality, the polymer molecules are stretching during the adsorption pro- 
cess, creating a polymer brush (for semitelechelic polymers) after sufficiently long time 
[TTJ |3] . Thus each molecule "covers" a smaller surface area at later times. Consequently, 
it is possible to adsorb more molecules onto the surface. Here, we take this fact into 
account and we modify the random sequential adsorption algorithm accordingly. We 
state our generalized random sequential algorithm (gRSA) in one dimension as follows. 

gRS A algorithm: We consider adsorption of small intervals of different sizes 
onto the interval [0, 1], the adsorbing domain. At each time step, the size of a 
small interval is chosen randomly according to the probability distribution func- 
tion p(z). We select randomly the position of the center w inside the adsorbing 
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domain [0, 1] and we make an attempt to place the small interval of length z 
there. If the center w of the segment to be adsorbed lies inside a segment al- 
ready placed, the adsorption is rejected. If the position of the center w is chosen 
in the gap (xi,x 2 ) between attached polymers, then the new polymer segment 
is adsorbed with probability £(z, w — Xi, x 2 — w) where w — X\ and x 2 — w are 
distances of the center w of the polymer from the endpoints of the gap (x\,x 2 ). 

The parameters of gRSA which have to be specified include the probability distribution 
function p(z) and the probability w — x%,x 2 — w). In what follows, we assume that 
the lengths of polymers are distributed according to the formula 



p(z) = 



+ l)e~ a - l z a 




for z < e, 
for z > e, 



(3.1) 



for a > — 1 and small e<l. Let us assume that the position of the center w of the new 
polymer is chosen in the gap [xi,^], i.e. w G [xi,^]- Then we take the probability 
(per unit time) of adsorbing the polymer segment of length z < x = X\ — x 2 as 
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and the probability of adsorbing the polymer segment of length z > x as 
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(3.3) 



In the latter case, the maximum probability of adsorption is achieved for w = X1 + X2 ; 
for which £(w) = -. The formulas (|3.2j) and ()3.3|) give the same probability density 
function ^(-) for z = x as is desirable. The plot of ^ as a function of w is given in Figure 
Formula (|3.2|) is shown in Figure C^a) where the gap size x = x 2 — x± is greater than 
the length of the new polymer segment z. Formula (J3.3)) is shown in Figure C^b) where 
the gap size x = x 2 — x% is less than the length of the new polymer z. 

To explain the motivation behind formula ()3.2|) . three possible cases of the relative 
position of the new (red) interval of the length z < x and the gap (xi,x 2 ) are shown in 
Figure|21 In Figure|2^a), the red interval is rejected because its middle point w lies inside 
a polymer segment which is already adsorbed to the surface. Hence, the probability of 
adsorption is 0, the same probability as in the classical RSA model. In Figure Wi c ), 
the red segment of the length z does not overlap with neighbouring polymers, and we 
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Xj x 2 w x 1 x 2 w 

Figure 1: The probability £ as a function of w for gRSA model (jj.jjj) - (jff.ffi : (a) for the 

case z < x = X2 — x±; (b) for the case z > x = X2 — Xi- 




w w w 

Figure 2: Schematic of gRSA. (a) Polymer is refused; (b) polymer is adsorbed with the 
probability £(z,w — Xi, x 2 — w); (c) polymer is adsorbed. 



allow it to be adsorbed with probability £ = 1. Cases in Figure E(a) and Figure 12(c) are 
treated as in the classical RSA model. 

In Figure 12(b), the center of the red polymer is inside the gap but the red polymer 
overlaps with neighboring polymers segments. This polymer would be rejected by the 
classical RSA model. We believe it is more realistic to consider that such a polymer 
will be adsorbed with some nonzero probability which continuously interpolates between 
the cases shown in Figures |2( a) and [2(c), i.e. between zero for w — x\ =0 and 1 for 
w — X\ = z. Formula (|3.2j) takes this fact into account, using simple linear interpolation. 
Formula ()3.3|) naturally extends the formula (|3.2j) for polymer segment lengths greater 
than the gaps (see also Figure HJ). 

Having explained the new rules for adsorbing the polymer, we must also specify what 
part of the surface is actually covered. We will assume that the new polymer covers 
only the intersection of the intervals 

This guarantees that a possibly long, newly adsorbed polymer will not "spill over" and 
cover any part of the neighboring gaps. 
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4 Analysis of gRSA 

Let G(x,t) be the concentration of gaps (holes) of length x at time t and let C(x,t) be 
the corresponding cumulative probability distribution function; that is, 



C{x,t) 



I™G{y,t)dy 



G(y,t)dy. 



(4-i; 



The total length of the surface that is covered by polymers at time t, A(t), is directly 
related to G(x, t) by 

A{t) = 1 - / xG(x,t)dx. (4.2) 
Jo 

The number of polymers attached to the surface at time t, N(t), can be also related 
to G(x,t), as we will see in Section T4.ll Thus, the starting point of the analysis of the 
system is the derivation of the evolution equation for the distribution function of gaps 
G(x,t). 

A gap of length x can be created from a larger gap (of length y > x) by adsorbing a 
suitable interval to the system. Thus the evolution of the concentration of gaps G(x,t) 
is given by the equation 



dG 
~dt 



(x,t) 



G{x,t) 



+ 





2(y-x) 



£(z, u, x — u)du 



p(z)dz + 



(4.3) 



2 £ (z, x + ~, y - x - |j p(z)dz 



G(y,t)dy. 



Using ()3.2|) and ()3.3|) . equation ()4.3|) can be rewritten in the following form 
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2{y-x) 



p{z)dz+ (4.5) 



G(y,t)p(z)dzdy. 



We assume that the lengths of polymers are distributed according to formula (j.llj) for 
a > — 1 and small £<1. Moreover, we assume that there are already no holes of the 
length greater than e/2 in the system, i.e. G(x,t) = for x > e/2. Then (using (jSHJ)), 
equation ()4.5|) can be rewritten (for x < e/2 and a ^ 0) as 



dG G(x,t)(a + 1 ) f x 

~dt [X,t) ~ i^+i 



x — 
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y-x 
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which implies 



dG x a+2 G(x,t) x 2 G(x,t)(a + l) 2 a+2 - 4 

dt ^' ' ~ a{a + 2)e a+l 2^e + ae a+l 



If a = 0, equation (|4.5jl implies (for x < e/2) 
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z a dzdy 



G(y,t)(y -x) a+ dy. (4.6) 



G(y,t)( 2 /-x)d 2/ . (4.7) 



Equation (|4.6|) (or ()4.7jl ) is the desired integro-differential equation for G(x,t). If we 
solve (|4.6|) . we can compute the evolution of A(£) by ()4.2|) . The equation for the evolution 
of N(t) is given in the next section. 



4.1 Evolution of N(t) 

At each time step, an interval of length between (z, z + dz) is chosen with probability 
p(z)dz. This interval can be placed in any gap of size x with probability J* £(z,u,x — 
u)du. There exist G(x,t)dx gaps whose size lies in the interval (x,x + dx). Hence, the 
integral J^°[Jq £,{z, u,x — u)du]G(x, t)dx gives the probability that the randomly chosen 
position of the polymer of length z will be accepted. Thus the probability of attaching 
a polymer of any length at one time step is equal to 



o L^o 



£(z, u, x — u)du 



G(x, t)p(z)dxdz. 



Using a continuous approximation for N(t), we find that N(t) satisfies the following 
ordinary differential equation 



dN 
~dt 



oo ^oo 



JO 
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(4.9) 
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Taking p(z) to be given by (jH.lj) and (,(z,u, x — u) to be given by (jH.2j) - (jH.Hj) . and 
considering the regime where all gaps are already less than e (i.e. G(x, t) = for x > e), 
we obtain 
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Hence 



dN 
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G(x, t)x 2 dx. 



(4.10) 



Before analyzing (|4.6|) and (J4.1U)) further, we summarize some results from the literature 
on classical RSA. 



4.2 Some results for the classical RSA 

If we consider particles of the same length e so that p(z) = 5{z — e), and if we choose 
£(z, w — xi, x 2 — w) of the form 

{z z 
1 for x 1 + - < w < x 2 - -; 
2 2 (4.11) 

otherwise; 

then our gRSA algorithm reduces to the classical RSA algorithm. The evolution equation 
(14. 3|) can be used to verify known one-dimensional results about fixed segment size, non- 
overlapping random sequential adsorption JUj, namely one can show that the jamming 
limit is asymptotically approached as t~ l p~3] . 

Random sequential adsorption with a probability distribution p(z) given by (jH.lj) and 
probability w — xi, x 2 — w) given by (j4.11j) has been studied in [TOj. Then equation 
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()4.3j) for x < e reads as follows (assuming that initially there exist no holes of length 
greater than e in the system, i.e. G(x, t) = for x > e) 

dG. , x a+2 G(x,t) 2 [ x+ c 

= - 7Z I n +-^ZT / G(y,t)(y-x) a+1 dy. (4.12) 



dt v ' y (a + 2) £ <*+ 1 
The scaling ansatz ^0] for the concentration G(x,t) can be written as 

G(x,t) ~ t a $ (xt 6 ) for £>1, and xt b finite. (4.13) 

Defining the moments 

M 7 (t) = / x 7 G(x,t)dx, m 7 = / £ 7 $(0 d £ ( 4 - 14 ) 
and using ()4.13j) . we obtain 

M 7 (t) ~ t a ^ 7 m 7 . (4.15) 

Moreover, multiplying equation (|4.12j) by x 1 and integrating over x, one can derive the 
equation for moments, 

^(^) = ^^ M 7+-+2 where F( 7 , a) = 25( 7 + 1, a + 2) - — ^ (4.16) 

where •) is Beta function. We define the function 7(a) implicitly by the equation 
F(j, a) = 0. If 7 is equal to 7(a), then the moment M 7 is independent of time. Hence, 
using (|4.15jl . we obtain the relation a = 6 + 67(0;) between the coefficients of the scaling 
ansatz ()4.13|) and the parameter a of the model. Finally, substituting the scaling ansatz 
(14.13)) in (|4.12|) . we find that b = (a + 2) _1 . Thus, the scaling of moments (|4.15|) can be 
rewritten in the form 

M p (t)~ir where /i — ( 4 - 17 ) 
Using (JOI) and (j4~T7|) . we obtain 

\-A(t)= [ xG(x,t)dx ~t-" {a) where uo(a) = 1 - ^ . (4.18) 
io « + 2 

The graph of the function u(a) is given in FigureOfa). The equation ()4.9|) for p(z) given 
by (|3.1j) and probability £(z, w — Xi, x% — w) given by ()4.11|) reads as follows: 

— = — - / x" +2 G(x,t)dx. (4.19) 
dt a + 2 J 

Using ()4.17|) . we obtain (for a (a) > 0) 

N(t) ~ where tr(a) = ^1. (4.20) 

The graph of the function a (a) is given in Figure Ofb). 
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Figure 3: (a) The graph of the exponent oj(a) given by l \4-18\j . (b) The graph of the 
exponent a(a) given by <\4-2CQ . 

4.3 Evolution of gRSA 

The temporal evolution of the gRSA model is more complex than the cases discussed 
in Section T4.21 To see this, we use the moments M 1 (t) defined by (j4.14j) . Multiplying 
equation 1)4.6)1 by x 7 and integrating over x, we can derive the equation for moments 
(for 

. . 1 xr a + 1 „ , 

coo oq+2 a roc 



a(a 



rco na+2 A roc 

+ jf a£a+1 J G(y,t)x^y-xr +1 dydx = 
1 cy 4- 1 Z" 00 2 Q+2 — 4 /" y 



-^7+a+2 



a(a + 2)e 



1 /2 a+2 — 4 1 \ a + 1 

-B( 7 + 1, a + 2) + — oN ) M 7+Q+2 - ^-^M 7+2 



V a w ' y afa + 2) J r ™ 2ae 

r #( 7 , «)M 7+Q+2 - ^M 7+2 , (4.21) 



v " ' 2a£ 
where £(•, •) is Beta function and #(7, a) is defined as 



( 7 , a) = 2 " +2 4 i?( 7 + 1, a + 2) + 1 (4.22) 
a a (a + z) 
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First, consider the case a < 0; the dominant term on the right-hand side of 1)4.21)1 is the 
term e~ a ~ 1 H( / -f, a)M 1+a+ 2. The zeroth order moment, 

poo 

M (t) = / G(x,t)dx, 
Jo 

gives the total number of gaps at time t. At leading order, we have 

^( x ,t)=e- a - 1 H(0,a)M Q+2 . 

There is a constant a = —0.415 such that H(0, a) is positive for a < a and negative for 
a > a. We immediately see that 

G(x,t)dx ^ for a > a. (4.23) 

To illustrate the result (|4.23jl . we will execute two stochastic simulations with the gRSA 
algorithm. We will use 1)3.1)1 . ()3.2)) and ()3.3)) where a = —0.1 or a = —0.3. We choose 
e = 10~ 3 . The results are given in Figure where the time evolution of the number of 
gaps and the number of adsorbed polymers are shown. Note that we use a logarithmic 
scale on the time axis because the long-term dynamics are very slow. For a = —0.1, 
the stochastic simulation was stopped when 99.9999% of the surface was covered. For 
a = —0.3, the stochastic simulation was stopped when 99.9994% of the surface was 
covered. 

Next, we will study the behaviour of the system for a < a. Here, we will assume 
the scaling ansatz ()4.13jl . Differentiating 1)4.22)1 with respect of 7, we obtain 

dH 2 a+2 - 4 r 1 

— (7, a) = 5( 7 + 1, a + 2) ^0(7 + 1) - ^0(7 + a + 3) (4.24) 

£77 a L J 

where ipo is the polygamma function. For each a < a, the equation 

H(%a) = (4.25) 

defines implicitly the function 7(a). If 7 is equal to 7(a), then the moment M 7 is 
independent of time. Hence, using (14.15)1 . we obtain the relation 

a = b + 67(a) 

between the coefficients of the scaling ansatz ()4.13)1 and the parameter a of the model. 
Finally, substituting the scaling ansatz ()4.13)1 into ()4.21j) . one can find that b = (« + 2) _1 . 
Hence, the scaling of moments ()4.15jl can be rewritten in the form 

Mp{t)~tP where JJ = (4-26) 
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Figure 4: gRSA model for a = —0.1 (panels on the left) and a = —0.3 (panels on the 
right). We plot the time evolution of the number of adsorbed polymers (top panels) and 
the time evolution of the number of gaps (middle panels). The time axis of the top and 
middle panels is logarithmic. We also plot the time evolution of the quantity [l — A^t)]^ 1 
(bottom panels) where time is scaled according to (4-27). 
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Figure 5: (a) The graph of the exponent U(a) given by !\4-27 ). The dashed line shows 
the exponent u (a) given by Q^.l<$p . (b) The graph of the exponent a(a) given by 
The dashed line shows the exponent a(a) given by ■ 



Using (JO]) and (j4^Hjl . we obtain 

1-A(t)= ( xG(x,t)dx~t~^ where 57(a) = 1 - ^ . (4.27) 
Jo a + 2 

The graph of the function u(a) is given in Figure Efa). We also plot u(a) given by 
()4.18|) for comparison. Using f)4. lUj) and ()4.26j) . we also obtain 

N(t) ~ f {a) where a (a) = (4.28) 

The graph of the function a (a) is given in Figure Efb); we also plot a(a) given by ()4.20|) 
for comparison. 

To illustrate the results (|4.27|) and (|4.28|) . we will execute two gRSA stochastic 
simulations. We will use (|3.1j) . (|3.2|) and (|3.3p where a = —0.5 or a = —2/3. We choose 
e = 10~ 3 . The results are given in Figure |H1 where the time evolution of the number of 
gaps and the number of adsorbed polymers are shown. The time is scaled according to 
()4.27|) and ()4.28|) ; we solve ()4.25|) to obtain the desired exponents 

a(-0.5) = 0.0872, u7(-0.5) = 0.5795, a ^-0 = 0.2875, W = 0.4625, 

(4.29) 

and then we plot the results of stochastic simulations using the corresponding scaling 
()4.29|) . In Figure |Hl we also plot the cumulative distribution function C(x, t) for 
different times (i.e. for different numbers of polymers attached). Using the suitable 
rescaling C(x,t) = C(kx,t), the curves collapse to a single curve as shown in Figure [7| 
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Figure 6: gRSA model for a = —0.5 (panels on the left) and a = —2/3 (panels on the 
right). We plot the time evolution of the number of adsorbed polymers (top panels) and 
the time evolution of the inverse gap size [1 — A(t)] -1 (middle panels). Time is scaled 
according to fl^.iffi ) (top and middle panels). We also plot the cumulative distribution 
function C(x,t) at different times for both simulations (bottom panels). 
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Figure 7: Scaled cumulative distribution function of gaps for gRSA model for a = —0.5 
(panel on the left) and a = —2/3 (panel on the right). 



Finally, one can easily show that formula (|4.27|) works also for the case a < a < 0. 
To illustrate it, we plot the time evolution of the quantity [1 — A(t)] -1 for a = —0.1 
and a = —0.3 in Figure 0] (bottom panel) using the scaling (|4.27|) . On the other hand, 
formula (J4.28)) is no longer correct because integrating of (|4.10|) implies that N(t) is 
approaching a constant value. More precisely, N(t) ~ C + t a ( a > for a (a) < 0. If a is 
positive than ()4.23|) is valid, i.e. we have ()4.23|) for any a > a. 



4.4 Time dependent concentration of reactive polymers 

As discussed before, the reactive groups on our polymers are capable of reacting with the 
solvent before reaching the surface jT^j. It may therefore be more realistic to consider 
that only a fraction of polymers r(t) is still reactive at time t. Here, r(t) € [0, 1], r(0) = 1 
and r(t) is a decreasing function of time. 

The random sequential adsorption algorithm has to be modified as follows: at each 
time step, we generate the random number uniformly distributed in the interval (0, 1). 
If this number is greater than r(t), then the selected polymer has lost its binding site 
through reaction with the solvent (it cannot be adsorbed), and we continue with the 
next step. Otherwise, we choose randomly a position on the interval and we attempt to 
place the polymer there. 

Depending on the form of function r(t), different dynamics can be observed. First, 
let us suppose that 

r(t) = \ for AG [0,1). (4.30) 
In this case, we can find a relation between the modified random sequential algorithm 
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and the previous results. At each time t, we can compute the average waiting time At 
before a reactive polymer hits the surface as the solution of the equation 



t+At -y 



-rdT = l. (4.31) 



Solving flOTJ), we find 

(t + Aty- x = t x - x + 1 - A. 



Hence, 

1/(1-*) 

At ' ' - ' 



+ 1 - A 



t x . 



Consequently, we can make use of the formulas ()4.18|) and ()4.2(Jj) . or formulas ()4.27|) 
and (|4.28jl . in the case (|4.30j) . For example, using (|4.18jl and (14.20)1 . we obtain that the 
quantities A(t) and N(t) satisfy the following asymptotic behaviour 

A{t) ~ t(!- A M«) and N{t) ~ t^- x) ^ a) . (4.32) 

To illustrate the formula 1)4.32)1 . we stochastically simulate the gRSA model with the 
probability distribution p(z) given by 1)3.1)) and the probability £(z, w — x%, x-i — w) given 
by ()4.11j) . where the fraction of the reactive polymers in the system decreases with time 
according to ()4.30j) . We select e = 10~ 3 and verify the asymptotic behaviour ()4.32)) for 
a = —0.5 and A = 0.5. Then ()4.32j) implies 

1 - A(t) ~ t-°' 1014 and N(t) ~ t 2319 . 

The time evolution of A(t) and N(t) is given in Figure |H1 (top panels). We also present 
results for a = and A = 0.33 in Figure |H1 (bottom panels). Again, we scale the time 
according to (J4.32)) . 

The reactive group is lost by chemical reaction with the solvent. It might be more 
natural to consider (instead of (J4.30)) ) that the fraction of reactive polymers is exponen- 
tially decreasing, i.e. 

r(t) = e~ xt for A > 0. (4.33) 

The formula ()4.33)) gives rise to qualitatively different dynamics for the system, as 
opposed to the dynamics associated with ()4.30j) . For simplicity, let us assume that every 
polymer with a functional reactive group can be adsorbed (which will give a bound on 
N(t) from above). Then the average number of adsorbed polymers N(t) is 



iV(t) = i(l- e - A ') ; 



which implies that the number of adsorbed polymers does not approach infinity as in 
the previous case. 
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Figure 8: Modified RSA model from Section \4-4\ The time evolution of the number of 
polymer molecules attached to the surface N(t) (left). Time evolution of the inverse of 
the total gap size [1 — Ait)]^ 1 (right). Time is scaled according to fl^.ffijp . 



5 Equation-free analysis of gRSA 

In the previous theory, we assumed the scaling ansatz ()4.13|) (see also [TU] ) for G(x,t) 
and we computed the time dependence of the quantities of interest A(t) and N(t). A 
related interesting question is whether we can also compute the profile $ from (|4.13|) . 
One possibility is to substitute ()4.13|) in equation ()4.6|) and solve it numerically for $ but 
we will not proceed this way. Instead, we demonstrate the computation of self-similar 
profile $ using only short-time appropriately initialized simulations of the stochastic 
gRSA model. In this equation-free context 012], it is easier to work with the cumulative 
distribution function C(x, t), which can be obtained from G(x,t) through (j4.1|) ; C(x,t) 
is less noisy than G(x,t) (e.g. jH]). Using (j4.1J) and (j4.13J) . we obtain 

1 r 

C(x,t) = „ / G(y,t)dy = 

Jo G{y,t)dy Jo 
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Hence, we see that the cumulative density function C(x,t) scales as 

C{x,t) = C(xt b ). 



(5.1] 



To compute the profile C, we can use an equation-free iterative fixed point algorithm 
which is shown schematically in FigureEl Starting with the initial guess Co, we compute 
the sequence of profiles Ck, K — 1,2,3, ... , where 

C K+1 = *(C K ), for K — 0,1,2,3, ... , (5.2) 

and where the mapping \1/ is obtained as the composition of the following four steps: 

(a) Given the cumulative density profile Ck, create one or more microscopic real- 
izations of gaps in the unit interval such that the initial cumulative density function 

is C(;0) = C K . 

(b) Use the microscopic simulator (i.e. use the gRSA algorithm) for a short time 
At. 

(c) Compute the new cumulative distribution function C(-,At) at time At. 

(d) Rescale C(-,At) to compute Ck+i- 

One possible way to rescale C{-, At) is to compute the average gap size a from C(-, 0) 
and the average gap size a& t from C(-, At). Then the Ck+i can be computed by 



C K+1 (x) = C[^x,At). 
a 



(5.3) 



We now present illustrative results obtained by this fixed point computation (|5.2|) using 
the gRSA algorithm. We will use (jSHJ), (Q and where a = -0.5 or a = -2/3. 
We choose e = 10~ 3 . The results of long term simulations for these parameter values 
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were already shown in Figure E3 Our goal is to use the iterative formula (|5.2jl to compute 
the scaled cumulative distribution function profile which was shown in Figure This 
algorithm allows us to find the self-similar shape by performing simulations while simu- 
lating at a scale (at relatively larger average gap sizes) where the evolution is relatively 
fast, compared to the long-term dynamics close to jamming. The initial guess is given 
as 

T , , _ j for x < 1.5 x 10" 4 ; 
o[X) ~\1 for x > 1.5 x 10- 4 ; 

which means that initially all our gaps have the same size 1.5 x 10~ 4 . At each iteration 
step (see Figure EJ), we place 1000 gaps according to the cumulative distribution function 
Ck to the interval [0, 1]. We evolve the simulation until 100 new polymers are placed. 
We then rescale the new cumulative distribution according to (|5.3jl and we compute 
Ck+i- Several first iterations are shown in Figure [TU] (top panels). We see that after 
20 iterations, we have effectively reached the steady state (the stationary shape of the 
self- similarly evolving gap distribution). More precisely, the error between iterations is 
small and it is not further systematically decreasing. The comparison of the equation- 
free 20th iteration with the results obtained by the long-time simulations are also shown 
in Figure ITU1 (bottom panels). 

Finally, we note that many other algorithmic possibilities for the computation of the 
profile C exist. The equation (|5.2|) seeks a fixed point of the mapping Instead of suc- 
cessive substitution, other fixed point algorithms implemented in a matrix-free fashion 
through short simulation bursts can be used to find stationary solutions - for example, 
Newton-GMRES iterations 2J; these would be able to converge on even dynamically 
unstable self-similarly evolving distributions. 



6 Discussion 

In this paper, motivated by a pharmacological example involving polymer coating of a 
virus surface, we studied certain aspects of polydisperse adsorption of macromolecules 
in one spatial dimension. We presented an extension of the classical random sequential 
adsorption algorithm to capture better certain essential properties of the pharmacolog- 
ical model system currently used in drug development research. We introduced partial 
overlapping of adsorbing macromolecules, i.e. we considered that the polymers are not 
rigid objects but they can be deformed while attaching to the surface. We found two 
distinct asymptotic regimes. Depending on the parameters of the processes involved, 
we can observe that either (a) the number of gaps between polymers asymptotically 
approaches zero, or that (b) the number of gaps asymptotes to infinity and the gap 
distribution acquires an asymptotically self-similar profile. 

We also briefly discussed the impact of a possible reaction of the polymers with 
the solvent on gRSA dynamics. Again, two possibilities exist. If the decay of the 
reactive groups is relatively weak, then the dynamics of the system remains qualitatively 
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Figure 10: Equation free gRSA computational results for a = —0.5 (left panels) and a = 
—2/3 (right panels). Iterations of the equation-free dynamic renormalization algorithm 
(top panels). Comparison with the steady steate profile obtained through long time 
simulations (bottom panels). 
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unchanged and the system only evolves on a slower time scale. On the other hand, if 
the reactive groups decay exponentially, this decay ultimately wins over the polynomial 
time asymptotics of gRSA. From the applications point of view, it therefore becomes 
crucial to know the corresponding rate constants in order to reliably predict what type 
of behaviour one might expect over the time scales of interest. Typically, the coating 
process is performed overnight in the laboratory and different reactive groups have 
different half lives; measuring these rates becomes an important task. 

In this paper we worked in one spatial dimension and provided analytical results 
about the long time behaviour of gRSA models. The analytical approach was based on 
two important facts: we knew what the good macroscopic observables for describing the 
system behaviour were and we were able to write down analytically tractable equations 
for these observables. The good observable for our system was a distribution of gaps 
G(x,t) between adsorbed polymers. If we know the initial distribution of gaps G(Q,t) 
one could easily predict G(x, t) at future times. 

On the other hand, if we know (hope) that the gap distribution G(x, t) is a good 
observable for the system of interest but we do not know the evolution equation for 
G(x,t), then it is still possible to use the equation- free methods The main idea 
of the equation-free methods is to use the short bursts of appropriately initialized mi- 
croscopic/stochastic computations to estimate macroscopic quantities of interest on de- 
mand. Hence, if one does not have an explicit coarse-grained evolution equation for the 
system statistics, one can in principle avoid long, brute-force simulations. This might be 
the case for one-dimensional adsorption problems with more complicated microscopic 
evolution rules. 

The situation becomes significantly more difficult in the higher dimensional case. 
Here, the analytical theory is far behind in development, and the literature contains 
mostly computational results. The first question for higher dimensional adsorption 
is the nature of the "good" coarse-grained observables for the system. Good observ- 
ables (the variables in terms of which the unavailable effective model would be written) 
are necessary for developing a useful analytical theory. Knowing appropriate coarse- 
grained observables is also an important feature of equation-free algorithms. Having 
one-dimensional analogues in mind, we see that one needs an effective way to describe the 
statistics of "gaps" (free space) in higher dimensions. If we cannot estimate (by intuition 
or by suitable algorithms for the detection of low-dimensionality in high-dimensional 
data) effectively good observables for the system, then the direct, brute-force compu- 
tationally intensive simulations might be the only modelling option. In this paper we 
showed cases where we could do more than brute-force simulation and provided analyt- 
ical results giving insights into the dynamics of gRSA. 

The problems studied in this paper were motivated by the pharmacological exam- 
ple mentioned above, and realistic predictive modelling of the problem clearly requires 
extensive model parameter information that must come from experimental data. As 
we showed, we can expect different dynamics of the problem depending on the values 
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of the parameters of the polymer and the virus which are used. Obtaining reliably 
such parameters and bounds on their uncertainty for our particular model problem is 
non-trivial, and we are not yet ready to report about it. 
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